% 三层神经网络反向传播示例 
clear; clc;

%% 1. 初始化参数 - 修正维度
% 网络结构：输入层(1) → 隐藏层(2) → 输出层(1)
W1 = [0.1; 0.2];        % 隐藏层权重 (2×1)
b1 = [0.3; 0.4];        % 隐藏层偏置 (2×1) - 修正为列向量
W2 = [0.5, 0.6];        % 输出层权重 (1×2)
b2 = 0.7;               % 输出层偏置 (标量)
x = 1;                  % 输入 (标量)
y = 1;                  % 真实标签 (标量)

fprintf('网络参数:\n');
fprintf('W1 = [%.1f; %.1f]\n', W1(1), W1(2));
fprintf('b1 = [%.1f; %.1f]\n', b1(1), b1(2));
fprintf('W2 = [%.1f, %.1f]\n', W2(1), W2(2));
fprintf('b2 = %.1f\n', b2);
fprintf('x = %.1f, y = %.1f\n\n', x, y);

%% 2. Sigmoid函数及其导数
sigmoid = @(z) 1./(1 + exp(-z));
sigmoid_der = @(a) a .* (1 - a);  % 更稳定的计算方式

%% 3. 前向传播
z1 = W1 * x + b1;       % (2×1) = (2×1)×(1) + (2×1)
h = sigmoid(z1);        % 隐藏层激活 (2×1)
z2 = W2 * h + b2;       % (1×1) = (1×2)×(2×1) + (1×1)
y_hat = sigmoid(z2);    % 输出层激活 (标量)
L = 0.5 * (y - y_hat)^2; % 损失 (标量)

fprintf('=== 前向传播结果 ===\n');
fprintf('z1 = W1*x + b1 = [%.1f; %.1f]*%.1f + [%.1f; %.1f] = [%.4f; %.4f]\n', ...
        W1(1), W1(2), x, b1(1), b1(2), z1(1), z1(2));
fprintf('h = sigmoid(z1) = [%.4f; %.4f]\n', h(1), h(2));
fprintf('z2 = W2*h + b2 = [%.1f, %.1f]*[%.4f; %.4f] + %.1f = %.4f\n', ...
        W2(1), W2(2), h(1), h(2), b2, z2);
fprintf('y_hat = sigmoid(z2) = %.4f\n', y_hat);
fprintf('损失 L = 0.5*(y - y_hat)^2 = 0.5*(%.1f - %.4f)^2 = %.6f\n\n', y, y_hat, L);

%% 4. 反向传播
% 输出层梯度
dy_hat_dz2 = sigmoid_der(y_hat);    % sigmoid导数 (标量)
dL_dy_hat = y_hat - y;              % 损失对输出的梯度 (标量)
dL_dz2 = dL_dy_hat * dy_hat_dz2;    % 损失对z2的梯度 (标量)

dL_dW2 = dL_dz2 * h';               % 损失对W2的梯度 (1×2)
dL_db2 = dL_dz2;                    % 损失对b2的梯度 (标量)

% 隐藏层梯度
dL_dh = W2' * dL_dz2;               % 损失对h的梯度 (2×1) = (2×1)×(1)
dh_dz1 = sigmoid_der(h);            % h对z1的梯度 (2×1)
dL_dz1 = dL_dh .* dh_dz1;           % 损失对z1的梯度 (2×1)

dL_dW1 = dL_dz1 * x;                % 损失对W1的梯度 (2×1) = (2×1)×(1)
dL_db1 = dL_dz1;                    % 损失对b1的梯度 (2×1)

%% 5. 输出结果
fprintf('=== 反向传播梯度计算过程 ===\n');
fprintf('输出层梯度:\n');
fprintf('  dL/dy_hat = y_hat - y = %.4f - %.1f = %.4f\n', y_hat, y, dL_dy_hat);
fprintf('  dy_hat/dz2 = y_hat*(1-y_hat) = %.4f*(1-%.4f) = %.4f\n', y_hat, y_hat, dy_hat_dz2);
fprintf('  dL/dz2 = dL/dy_hat * dy_hat/dz2 = %.4f * %.4f = %.4f\n', dL_dy_hat, dy_hat_dz2, dL_dz2);
fprintf('  dL/dW2 = dL/dz2 * h'' = %.4f * [%.4f, %.4f] = [%.4f, %.4f]\n', ...
        dL_dz2, h(1), h(2), dL_dW2(1), dL_dW2(2));
fprintf('  dL/db2 = dL/dz2 = %.4f\n\n', dL_db2);

fprintf('隐藏层梯度:\n');
fprintf('  dL/dh = W2'' * dL/dz2 = [%.1f; %.1f] * %.4f = [%.4f; %.4f]\n', ...
        W2(1), W2(2), dL_dz2, dL_dh(1), dL_dh(2));
fprintf('  dh/dz1 = h.*(1-h) = [%.4f; %.4f].*[%.4f; %.4f] = [%.4f; %.4f]\n', ...
        h(1), h(2), 1-h(1), 1-h(2), dh_dz1(1), dh_dz1(2));
fprintf('  dL/dz1 = dL/dh .* dh/dz1 = [%.4f; %.4f] .* [%.4f; %.4f] = [%.4f; %.4f]\n', ...
        dL_dh(1), dL_dh(2), dh_dz1(1), dh_dz1(2), dL_dz1(1), dL_dz1(2));
fprintf('  dL/dW1 = dL/dz1 * x = [%.4f; %.4f] * %.1f = [%.4f; %.4f]\n', ...
        dL_dz1(1), dL_dz1(2), x, dL_dW1(1), dL_dW1(2));
fprintf('  dL/db1 = dL/dz1 = [%.4f; %.4f]\n\n', dL_db1(1), dL_db1(2));

fprintf('=== 最终梯度结果 ===\n');
fprintf('dL/dW2 = [%.6f, %.6f]\n', dL_dW2(1), dL_dW2(2));
fprintf('dL/db2 = %.6f\n', dL_db2);
fprintf('dL/dW1 = [%.6f; %.6f]\n', dL_dW1(1), dL_dW1(2));
fprintf('dL/db1 = [%.6f; %.6f]\n', dL_db1(1), dL_db1(2));

%% 6. 梯度验证（可选：数值梯度检查）
fprintf('\n=== 梯度数值验证 ===\n');
epsilon = 1e-6;

% 对W2的数值梯度
numerical_dW2 = zeros(size(W2));
for i = 1:numel(W2)
    W2_plus = W2; W2_plus(i) = W2_plus(i) + epsilon;
    z2_plus = W2_plus * h + b2;
    y_hat_plus = sigmoid(z2_plus);
    L_plus = 0.5 * (y - y_hat_plus)^2;
    
    W2_minus = W2; W2_minus(i) = W2_minus(i) - epsilon;
    z2_minus = W2_minus * h + b2;
    y_hat_minus = sigmoid(z2_minus);
    L_minus = 0.5 * (y - y_hat_minus)^2;
    
    numerical_dW2(i) = (L_plus - L_minus) / (2 * epsilon);
end
fprintf('数值梯度 dL/dW2 = [%.6f, %.6f]\n', numerical_dW2(1), numerical_dW2(2));
fprintf('解析梯度 dL/dW2 = [%.6f, %.6f]\n', dL_dW2(1), dL_dW2(2));